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Phase correlations, density fluctuations and three-body loss rates are relevant for many 
experiments in quasi one-dimensional geometries. Extended mean-field theory is used to 
evaluate correlation functions up to third order for a quasi one-dimensional trapped Bose 
gas at zero and finite temperature. At zero temperature and in the homogeneous limit, we 
also study the transition from the weakly correlated Gross-Pitaevskii regime to the strongly 
correlated Tonks-Girardeau regime analytically. We compare our results with the exact Lieb- 
Liniger solution for the homogeneous case and find good agreement up to the cross-over 
regime. 
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I. INTRODUCTION 

Low dimensional physics has always been in the focus of interest as reduced dimensionality 
goes hand in hand with an increase of quantum fluctuations. Thus, by reducing the volume of 
accessible phase space, the effect of quantum fluctuations of the remaining degrees of freedom 
must be enhanced in order to comply with the fundamental Heisenberg uncertainty principle. 

Currently, this fact has stimulated many fascinating experiments in the context of ultracold 
gases [[Tl|2l[3l|4l|5l[6l|71[8l, which explore various aspects of the geometric transitions. Serendip- 
itously, also most exactly solvable models of field theory are one-dimensional [9] and rest on 
the celebrated Bethe-ansatz invented in the 1930ies. In the context of atomic Bose gases, today 
most prominent are the spatially homogeneous models of hard- and soft-core bosons on a string 
of M. Girardeau [TOl as well as E. Lieb and W. Liniger [fTTl \T2\ . One of the many interesting 
questions which can be explored, is the cross-over from the weakly correlated Gross-Pitaevskii 
regime (7 <^ 1) to the strongly correlated Tonks-Girardeau regime lfT3l [T4ll (7 ^ 1). Thereby, 
one commonly uses the Lieb-Liniger (LL) parameter 7, cf. (|5]), to measure the relative strength of 
kinetic to repulsive self-energy in the dilute Bose gas. 

Today, we also have alternative tools available: First, there are the exact few-body calcula- 
tions, i. e. multi-channel time-dependent Hartree-Fock (MCTHF) or configuration interaction (CI) 
methods, which originate from atomic, molecular and nuclear physics. While originally designed 
for fermionic energy structure calculation, they are nowadays also applied to few-boson systems 

10-100 particles) in arbitrary trap geometries [fT5l [T6l [TTl [T8l . Second, there is now the pos- 
sibility to prepare atomic gases in optical lattices and this opens up the rich methodology of the 
Bose-Hubbard model and density matrix renormalization group methods lfT9ll20ll2ni22ll23ll24| . 
Third, there are stochastic multi-mode trajectory simulations [|25l that also successfully address 
the same questions. 

Irrespective of the choice of method, all need to predict experimentally accessible observables 
in terms of correlation functions, spatial averages or Fourier transforms, thereof. Most relevant are 
obviously the lowest order moments of the bosonic field operator a^, which is the single particle 
density rix at position x and the conjugate phase quadrature correlation function g^xl- The fluc- 
tuations about the mean density are measured with the second order density-density correlation 
function (^i^y 

/^t- \ (1) (^i^x) (2) (^x^l^j/^x) 
nx = {alaj, g)^l = ' , g^f^' = "-^ . (1) 
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In here, (...) = Tr{. . . p} denotes an average over the state of the system described by the many- 
body density operator p. Such second order correlation functions have been measured experimen- 
tally ||4l|51|26l|271|28l|29ll, while the third order density-density-density correlation 

(3) _ {^l^l^l^z%^x) 
9x,y,z ~ ' W 

became observable only recently [I3l|30l via the three-body recombination rate ||3T1l . 

Theoretically, much attention has been directed towards second order correlation functions 
Il32l[33ll34l[35ll36l[37l[38ll39l . while less is known about the third order correlation function. This 
situation has been rectified recently in [|40l |4T1| where the diagonal behaviour of this correlation 
function was calculated in the framework of Lieb-Liniger theory. This is where extended mean- 
field theory is useful, because we can calculate arbitrary orders of the correlation function and 
will present calculations of the diagonal and off-diagonal behaviour of the third order correlation 
function at zero as well as finite temperature. However, the extended mean-field (EMF) approach 
is restricted to values of the correlation parameter 7 < 1, because any mean-field theory is known 
to fail in the strongly correlated regime. 

This paper is organized as follows: In section |ll| we briefly review the central ideas of Lieb- 
Liniger theory [[TTll . This celebrated solution of the one-dimensional homogeneous Bose gas is 
an ideal benchmark for the extended mean-field theory |l36l |42l |43l |44l, whose basic concepts 



are summarized in section III In section IV we specialize the kinetic equations to a quasi one- 
dimensional homogeneous situation at zero temperature for which analytical solutions can be 
found and compare correlation functions with the LL-predictions. The extension to inhomoge- 
neous, harmonically trapped systems at finite temperatures is studied numerically in section [V} 
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II. LIEB-LINIGER THEORY FOR BOSONS IN ONE DIMENSION 

Lieb-Liniger theory based on the Bethe ansatz [|9l describes a one-dimensional homogeneous 
gas of N bosons on a ring of length L. It is one of very few exactly solvable problems in many- 
body physics and provides a solution for every value of the correlation parameter 7. Even in 
inhomogeneous trapped systems this is very useful, if we can make the local density approxi- 
mation. In the language of second quantization, the starting point for Lieb-Liniger theory is the 
following Hamiltonian 

H = dxall - — TT^ + -ala,^ a^, (3) 



Jo V 2m ' 2 ^ y ^' 

where m denotes the mass of a boson and the creation and annihilation operators satisfy the usual 
bosonic commutation relation. With the help of the Hellmann-Feynman theorem [|45l . one can 
obtain the diagonal part (x = y = 0) of the translation invariant second order correlation function 
gi^l = g^l, introduced in ([l]), as 

by differentiating the ground state energy Eq with respect to the coupling constant g. Here, |^o) 
represents the ground state and n = N/ L denotes the linear particle density. It was shown by 
Lieb and Liniger [1 IJ that the ground state energy only depends on the dimensionless correlation 
parameter 7. It is basically the ratio of the repulsive mean-field energy gn to the kinetic energy 
h^/2md'^ at an average distance d = 1/n. Another length scale of the problem is the healing 
length ^, which equates the kinetic energy of a wave function at scale ^ to the mean-field energy 

= — ^ = ^ (5) 

fi^n^ y/2'mng 

We call bosons weakly correlated for 7^1 (Gross-Pitaevskii regime) and strongly correlated for 
7^1 (Tonks-Girardeau regime). 

In terms of this parameter, the ground state energy and second order correlation function 



Eo = N'-^e{j), ^?S = e'(7), (6) 
are given in terms of the solutions of the Lieb-Liniger equations 

e(7) = -^^^ J dxh{x,'j)x (7) 
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FIG. 1: The Lieb-Liniger correlation function versus the correlation parameter 7. In subplot a), we depict 

(2) (2) 
the second order correlator g)^!^ (solid line), the Gross-Pitaevskii approximation g]^^ (dashed dotted 

line) and the Tonks-Girardeau approximation gj^^ j.^ (dotted line), while subplot b) shows the third order 

(3) (3) (3) 

correlator g)^!^ (solid line) and the approximations 5]^^ (dashed dotted line) and g)^!^ j.^ (dotted line). 

In the weakly correlated Gross-Pitaevskii limit 7 ^ 0, as well as in the strongly correlated 
Tonks-Girardeau regime 7 ^ 00, one obtains for the correlation function [1351 

&GP = 1 - ^> for 7 « 1, 9^lItg = ^> for 7 » 1. (9) 



A comparison of these approximations with the exact solution is presented in figure [T] The va- 
lidity of these results has recently been tested experimentally over a wide range of the correlation 
parameter [i5| and will be used to probe the extended mean-field approach presented in the next 
section. 
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FIG. 2: The auxiliary Lieb-Liniger function h{x, 7), as a function of the variable x for five different values 
of the correlation parameter 7. 

It is also possible to derive an exact result for the diagonal part of the third order correlation 
function within the framework of Lieb-Liniger theory, however the task is considerably more 
difficult. The exact result is derived in iHTTl by introducing a new function 6(7) which has the form 

.,5 rl 



e(7) 



7 



A^(7) 7-1 

and with the help of this function one obtains 



da; h{x, 



(10) 



(3) _ 36^7) - 4e(7) - 6e(7)eX7) ,r\ s , 9e^(7) - 5g(7) 

27 +\^+ ^2 



9ll 



(11) 



A comparison of the exact result for the third order correlation function with the approximations 
in the Gross-Pitaevskii and the Tonks-Girardeau regime [|35l is presented in figure [T] 



(3) 

9ll,gp 



^ for7«l, 9^lItg-Y^,: 



for 7 > 1. 



(12) 



The auxiliary function h{x, 7) of ^ is depicted for various values of the correlation parameter in 
figure |2] 
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III. EXTENDED MEAN-FIELD THEORY FOR BOSONS IN ONE DIMENSION 

A. Time-dependent Hartree-Fock-Bogoliubov equations 

The evolution of a weakly interacting dilute gas of bosons in three dimensions can be described 
by a Hamiltonian 



(13) 



p2 11 

Hxy = + Kxi(x)|y), Veccti^) = -mcu^x^ + -mculiy^ + z^), (14) 

where H is the single -particle energy in an external potential Vext and Hm is the two-particle 
potential. As we are interested in the quasi one-dimensional limit, we will consider a cigar shaped 
trapping configuration (angular frequencies to and uj±) with a large aspect ratio /?. The energy and 
length scales will be set by the transverse oscillator 

/3 = uj±/uj ^1, a_L = a/ h/muj±, e± = hu±. (15) 

Conceptually, it is straight forward in the extended mean-field theory to use real, finite range 
binary interaction potentials and obtain proper two-body T matrices including many-body correc- 
tions. However, for convenience, we will use the pseudo-potential approximation in here 

Hm(x - y) = dx-y, (16) 

m 

where denotes the s-wave scattering length. In order to compactify the notation we will inter- 
changeably use a subscript notation also for continuous functions. 

Extended mean-field theory uses a reduced state description based on a set of master variables 
{i E Basically, this implies the existence of a well separated hierarchy of time, energy and 

length scales [|46l |47| and leads to a rapid attenuation of correlation functions. Mathematically 
speaking, it allows for a selfconsistent expansion of the full many-body density matrix p in terms 
of a perturbation series of simple many-body density matrices cr^*), which depend parametric ally 
on the master variables 

p = a|°^) +a« (17) 

This non-perturbative series in terms of the interaction potential Vtin has been introduced first by 
Chapman and Enskog in the context of kinetic theory of gases [|48| . In addition to the simple series 
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expansion, we impose a selfconsistency constraint such that the operators 7^, corresponding to the 
c-number master variables 7j, fulfill 

7, = (7.) = Tr[7.p] = Tr[7,o-{°)j]. (18) 

As far as the master variables are concerned, we choose the mean-field a.^, the normal fluctua- 
tions of the single-particle density / and the fluctuations of the anomalous two-particle correlation 
function m such that 

(«x) = «x fi% = «x«y, = «x«y, (19) 

(OyOx) = /S + /x,y, (ay«x) = "^£y + ^K,y (20) 



A Gaussian operator cr|°| is compatible with the requirements of (|18 



19 



20). In turn, this im- 



plies the factorizability of multi operator products (Wick's theorem) and also yields non-Gaussian 
corrections by calculating the contribution of (^f^y 

By studying the coordinate transformation properties of the fluctuations [|49ll . one finds that the 
averages / and m are components of a positive semi-definite, generalized density matrix G > . 
Thus, the system is described by a row vector x, containing the mean-field a as well as its complex 
conjugate, and by the density matrix G 









Xx= 1 




, Gx.y - 1 



(21) 



/x,y '^x.y 
^x,y '^x,y + /x,y 

It can be shown from a Cauchy-Schwartz inequality that at T = 0, the generalized density matrix 
G obeys an idem-potency relation 

/i o\ 

GffsG + G = 0, (73 = , (22) 

Starting with the Heisenberg equation of motion, it is straightforward to derive the equations of 
motion for x and G. In order to obtain higher-order correlation functions within the present ap- 

2|), one has to evaluate the Gaussian Trf. . . crj^}] 



proximation scheme p2ll50l . like g^'^^ or g^^^ of ( 1 



as well as the non-Gaussian contributions Tr[. . . crj^}]- However it is clear that the Gaussian con- 
tribution will dominate for weak correlations. Thus, we have evaluated in here only the Gaussian 
contributions. But already for 7 1, deviations from that can be noticed. 
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B. Reduction to a quasi one-dimensional, stationary configuration 

In a very prolate trap, the transverse motion in the directions y and z is effectively frozen out 
and only amplitudes proportional to the ground state 

v>o{y,z,t) = — (23) 

Vvra_L 

need to be considered. By projecting all three-dimensional functions onto the longitudinal axis 
X ^ x, one obtains the time-dependent Hartree-Fock-Bogoliubov equations (THFB) for Xx and 

ihdtxx = lixXx + O(nL), imG.,.> = T.^G^^^, - h. c. + 0{VD, (24) 
with the following abbreviations for the single particle Hamiltonian and the self energies 





Ha 1 

























n. = , S,. = , (25) 



Hiv = + gf^ + 2(7/,,,, = + 2gj^ + 2g],^^, (26) 

In the course of the dimensional reduction, we had to introduce an effective one-dimensional 
coupling constant g = 2hiu±as, which is in agreement with previous derivations [|36l 15111521 . In 
order to obtain the stationary solution for the time-independent fields Xx and Gx,x', we make the 
ansatz 

X{t) = e-i^'^'^x, G{t) = e-s'^*'^3 ^gi^toa (28) 

which introduces the chemical potential fx and employs the Pauli matrix (J3 of ( 22 ). This ansatz im- 
plies that the normal fluctuations (t) become time-independent, while the anomalous fluctua- 
tions rhxy (t) oscillate with twice the chemical potential. The properties of the resulting stationary 
HFB equation will be investigated in section [IV| in the case of a homogeneous gas of bosons and 
in section |V] for a harmonic trapping potential, both at zero and for finite temperatures. 

The fact that the eigenvalue in the resulting stationary equations is indeed the chemical potential 
P9ll can be seen from a variation of the total energy functional £{a, /, m) = {H) given by 

£ = y dxdyS{x-y)[al{nx + '^\aJ^)a^ + {ny + gfxJfx,y] 

+ 1 j dx[(2|aj7.,. + + ^|m,,j2)+h.c.] + 0(\/l), (29) 

with the constraint that the number of particles N = J dx(/ifi + fx,x)- 
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IV. ANALYTIC SOLUTION FOR THE STATIONARY HFB-EQUATIONS IN THE 
HOMOGENEOUS SYSTEM AT ZERO TEMPERATURE 

For the calculations that are presented in the following sections, we use standard parameters 
for ^''Rb in natural units of length a± and energy e± 

m = 1.4432 ■ 10-25 kg, a, = 5.8209 ■ lO'^ m, 

uj± = 2n- 800 Hz, = 27r ■ 3 Hz, (30) 

a± = 3.8128 ■ 10-^ m, ^ = 2as/a± = 3.0533 ■ lO^^. 

To reach the homogeneous case, we have to decrease the influence of the external trapping 



potential in (27) by weakening it to the limit /3 ^ 1, where we can neglect it practically. There- 
fore, the equilibrium state should possess the same translation symmetry as the generator of the 
dynamics. Consequently, we can assume that the mean field is space independent = a and the 
density matrix only depends on relative differences r = x — x' 

/°° dk 
-7^e-'''''Qk. (31) 

Translation invariant systems are best described in Fourier-space, which was introduced above. We 
can also choose the mean-field to be real- valued by a suitable phase rotation. This is a consequence 
of the global number conservation that is built into the dynamical HFB equations [|53l . As the 
mean-field is real-valued f'^'^^ = m^^^ = , so are the fluctuations /o = fx,x and itiq = m^.^ and 
with these assumptions, the normalization constraint reads 

^=J = f^""^ + k (32) 

where n is the linear particle density on a length L. Furthermore the THFB equations ( [24] ) simplify 
significantly to 

/i = ^(/(^) + 2/o + mo), = (Sfc - ^a3)6?fe - h. c. (33) 

From the equation for the chemical potential, it is clear that energy and length scales emerge. 
It will be beneficial to introduce such scales for coherence {kc^c^^c), the pairing correlations 
(k, ^, Co), and their weighted sums and differences as 

K = = V^c, ~k = = v^, k^ = ^±'' = (34) 



u;, = A~gf^'\ u = -A~gmo, = 5^^^. (35) 
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In particular, in the Gross-Pitaevskii regime one can assume that ujc^ oj. With these definitions, 
one finds that the self energy is simply a 2 x 2 matrix in /c-space with eigenvalues ujk 




-k"^ — Co" I 



(36) 



Now, we can finally evaluate the density matrix part of the HFB equations (33 ). Moreover, we also 



have to consider the idem-potency relation of (22). It holds for the vacuum state at zero tempera 



ture and one obtains another, now quadratic relation between normal and anomalous fluctuations 
in A;-space 



rrik 



1 



fk = ml- fl 



(37) 



A;2 + CU+ V 2 

The system of equations can be solved point wise in A;- space and leads to two solutions. One of 
which has to be rejected on physical grounds. Thus, we find 

+ UJ+ 1 



rrik 



UJ- 



fk 



(38) 



The high momentum tail of the correlation functions is responsible for the short scale behaviour 
in real space. In the limit A; — > oo the leading terms are 



UJ- 



2k^ 



fk ~ ml. 



(39) 



A. Diagonal contributions of normal and anomalous fluctuations 



In order to obtain the diagonal part of the translation invariant correlation functions rhr and fr 



at r = 0, we have to evaluate the inverse Fourier transform of pi) 

cj„ r dk 



mo = -- 



u^k' 



fo 



dk~ 
27r^'" 



(40) 



Serendipitously, this can be done exactly in terms of elliptic integrals ||54 



fcc / ^ , mo\ mo 



/o = mo - 



27r 



(41) 
(42) 



where K and E are the complete elliptic integral of the first kind and second kind, respectively. 
Basic definitions are given in |A 1 
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FIG. 3: Diagonal part of the anomalous fluctuations —mo (left scale, solid line) and normal fluctuations fo 
(right scale, dashed line) versus mean field density f^'^K The asymptotic approximations for rfiQ (dashed 



dotted line) and /o (dotted line) according to ( 43 1 agree well for the considered parameter range 



The scaling properties of the correlation functions are most relevant for a physical insight. 
Thus, we can study the Gross-Pitaevskii regime of weakly correlated bosons, where | rrao//^^-* | <^ 1 
and use a series expansion for the elliptic integral (1 — x)K{l — x) ~ \n{4:/^/x). With this 
approximation we get 



UlQ 



W 64n J f(-)/~g , /, 



TT 



- rJiQ. 



(43) 



In this explicit formula, we had to introduce the Lambert- VT-function, which is defined in A 2 and 
an excellent asymptotic expansion is given in terms of logarithms [|55l . 

The approximations for the fluctuations are compared to exact numerical calculations in figure 
[3] and give a good agreement. We will use these approximations in the following sections to 
evaluate the ground state energy and correlation functions. 
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6. Off-diagonal contribution of normal and anomalous fluctuations 

1. Short length scale behaviour: r ^ ^ 

A rather simple, yet surprisingly efficient insight into the short range behaviour of the off- 
diagonal of the fluctuations can be obtained by using an iteration scheme for their Fourier trans- 



forms, which has its origin in (|37[). Starting with /^"^ = and using the recursion relation 



+ 



we get a rapid convergence towards the exact results. It is remarkable that even with the inverse 
Fourier transforms of low orders of this iteration scheme, we get a functional behaviour for /(r) 
and m(r), which is equivalent to their exact behaviour for short ranges. However in contrast to 



the exact form of the Fourier transforms of the fluctuations in (38), it is possible to perform the 
inverse Fourier transform analytically in every order of the iteration scheme. A closer look reveals 
that the dependence of the fluctuations on r has to be of the form 

m«(r) = e''=+l"lp«(r), /«(r) = e"'^+l"lQ(*)(r) (45) 

where P*^*)(r) and Q^'^\r) are polynomials in r of order 2' — 2 and 2*+^ — 3 respectively. Conse- 
quently the length scale on which the correlations decay is given by 

^+ = 7^ = ^ =. (46) 

k+ V2^(/W-mo) 

In the Gross-Pitaevskii regime, where f^"^^ ^ ttiq, we recover the healing length ,^ ^ 1/ a/2^/W, 
already introduced at the beginning in ([5]). 

The short range behaviour of the anomalous fluctuation and the normal fluctuation is depicted 
in figure |4j There, we compare the Ath order result of the iteration scheme to the exact numerical 
evaluation of the inverse Fourier transform. We assumed = 100 particles, distributed over a 
length of L = 90a_L. This length was chosen such that the density in the homogeneous case is 
similar to the density in the center of the trapped system, which will be discussed in section|Vj One 
obtains a good agreement between the approximation and the exact results in the regime where 
r <^ ^ ^ 10.5 with ^ ^ 3.88. At the origin we note that the anomalous fluctuation shows the 
typical cusp whereas the normal fluctuation has a smooth behaviour and consequently a vanishing 
first derivative at r = 0. 
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FIG. 4: Off-diagonal part of the normal and anomalous fluctuations versus the distance r. In subplot a) 
we depict the normal fluctuation /(r) (soUd line), its short range approximation f^^\r) (dashed dotted 
line) and the long range approximation (dotted line) according to ( [49] ). Subplot b) shows the anomalous 
fluctuation — m(r) (solid line), the short range approximation — m(^)(r) thereof (dashed dotted line) and 



the long range approximation (dotted line) according to (48 1. 



2. Long length scale behaviour: r S> ^ 



In order to get an approximation for the fluctuations in this regime, we start with the Fourier 



transform of the anomalous fluctuation in (38 1 and note for further consideration that the Fourier 



transform of the modified Bessel function of the second kind -ft'o(c|r|) is given by -n / \/ k'^ + . 
Thus the convolution property of the Fourier transform yields 



m[r] 



27r2 



KQ{k,\r'\)KQ(k\r' + r\)dr' . 



The modified Bessel function of the second kind -ft'o(^) is presented in figure [5} KQ{r) diverges 
logarithmically at the origin and decreases exponentially for large arguments 

where 7e ~ 0.5772 denotes Euler's constant. 



Koir) 





oo 



(47) 
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FIG. 5: The modified Bessel function of the second kind i^o(^) versus r. 

In the Gross-Pitaevskii regime where f^'^^ ^ l^^^ol^ we get from r ^ ^ that also r ^ and 
therefore the first Bessel function in the integral closely resembles a (5-function. Thus we obtain 
the result 



m{r) 
Rr) 



An 
Sirk 



1 + 



mo 



Ko [klr 



/(c) ^ 

^[klKoCk\r\)-k'K2Ck\r\ 



kc 
in 



Ko(k\r 



(48) 
(49) 



An alternative derivation of this result with the help of complex integration is given in|B] In figure 
|4]the asymptotic form of the fluctuations in terms of the Bessel functions is compared to an exact 
numerical simulation with the parameters that were mentioned in the previous subsection. In 
either case the semi-logarithmic plot reveals an exponential decay for large distances r i^, in 
agreement with the asymptotic behaviour of the Bessel functions. 

C. Comparison to Lieb-Liniger theory 



Having all the ingredients at hand to calculate correlation functions, we are now ready for a 
quantitative comparison with the results of Lieb-Liniger theory which provides an exact solution 
for the behaviour of the second and third order correlation function. 



As outlined in section III A we have to obtain the values of the correlation functions from 



multiple operator averages. If 6 is such a general operator then 



{7} 



(50) 



While we have already evaluated the Gaussian and non-Gaussian averages for the multinomial op- 
erator averages [|42l . it is clear that the Gaussian contribution will dominate for weak correlations. 
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Therefore we will focus in here on the Gaussian contribution and disregard the non-Gaussian 
contributions in the following explicit expressions of order two and one, respectively 



tlx,y 



.(2) 



.(3) 
JX,y 



y^flx fly 



+ 0{~g 



1 + 
1 + 
+ - 



y-'y,-- 



'^,y y,- 



n.r n. 



y 



'^^{fx,yfy,x ~^ ^x,y '^y,x) ~^ fx,yfy,x ~^ ^x,y'^y,x 



^x ^y 



+ 0{~g), 



(51) 

(52) 
(53) 



f(c) f If f I Of f 4- 2m* 

Jx,xJy,y ' J x,xJ y,y ' Jx,yJy,x ' x,y 



m. 



y,x 



+ - 



[m^ylmy^y) + rrly^yrhy^y + f y ^ f y ^y 



rir n. 



y 



x,yy^y,y^*y,x + ^H^lx) 



^^y 



^''^^^y,y^*y,x 



+ f f 



+ rniiiflym;,, + m;jyj + O(^), 



where rij. = fx'^l + f^^x denotes the total density. This way we can calculate the full diagonal 
and off-diagonal behaviour of the correlation functions. It works equally well for the trapped and 
homogeneous case. 

In figure|6]we see a comparison of approximations and exact numerical results within extended 
mean-field theory as well as Lieb-Liniger theory. In either case we observe a good agreement 
between our results and Lieb-Liniger theory. However as 7 increases the deviation from the exact 
result grows. We attribute this deviation to the non-Gaussian contributions that have been dropped. 

Another relevant quantity is the ground state energy of the system. By comparing the value of 



the energy functional ( 29 1 with the Lieb-Liniger ground state energy for a range of the correlation 
parameter 7, we obtain figure |7j In particular, we plot the relative deviation of the ground state 
energies. This is to be compared with deviations from a simple mean-field approach neglecting 
fluctuations and for the Bogoliubov method in the GP-regime which includes excitations of the 
mean-field. The latter approach results in [|56l 

4 

3^ 



e(7) 



LL,GP 



7 



3 

-72. 



(54) 



In either case, we present all the results in the form of a normalized deviation from the dimen- 
sionless ground state energy 6(7) from Lieb-Liniger theory given by (|6]). 




FIG. 6: Second and third order correlation functions versus the correlation parameter 7. In subplot a) we 
compare the exact results from Lieb-Liniger theory, gj^l (solid line), and the approximation in the GP- 



(2) (2) 
regime, (dashed dotted line) to gxA 



(2) 

q calculated with an extended mean-field theory. We 



depict exact results (dashed line) using (41 1 and (42) as well as approximated results (dotted line) using 



(43 1. In subplot b) we depict the same comparison for the third order correlation function. 




FIG. 7: Relative deviations from the dimensionless ground state energy of Lieb-Liniger theory as a function 
of 7. Results from an extended mean-field approach (solid line), the simple mean-field theory without 
fluctuations (dotted line) and the Bogoliubov approach (dashed dotted line) are compared. 
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The results show a clear improvement over simple mean-field theory and it also improves on 
the Bogoliubov method. Up to the cross-over at 7 1 the maximum deviation of our results is 
less than 4% and we obtain reliable results throughout the region of interest, i.e. 7 < 1. However, 
this appears to be the limit for a quasi one-dimensional extended mean-field theory and different 
approaches have to be used in the strongly correlated regime. 
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-100 



(c) 

FIG. 8: Coherent single particle density matrix fx,y versus x and y. As the ground state is real valued, the 

(c) 

coherent part of the pairing field rrix/y is also represented in this figure. 

V. NUMERICAL RESULTS FOR TRAPPED ATOMS AT ZERO AND FINITE TEMPERATURE 



A. The zero temperature limit for a trapped gas 



In the previous section we have studied the homogeneous case. In here, this will be extended 
to harmonically trapped systems and we present correlation functions up to third order. First of all 
we depict the spatial shape of the master variables /, m and of the quantities f^^\ m^^\ which are 
essential for the calculation of the correlation functions. The plots show numerical simulations for 



a particle number of = 10^ in a trap with standard parameters for ^''Rb according to (30). 

The coherent contribution to the single particle density matrix /i^ in figure 8 has off-diagonal 
long range order and extends over the complete system. As the Hamiltonian for a one-dimensional 
trap is real- valued, so is the ground state solution a^. Hence, the coherent contribution of the 
pairing field m^l is identical to fx^l and shown in figure [sj 

In contrast to the coherent contributions, the normal fluctuation f^^y in figure 9 and the anoma- 



lous fluctuation nix^y in figure 10 are primarily localized along the diagonal. 
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FIG. 9: Normal fluctuations f^^y versus x and y. 



The coherence in the off-diagonal direction is only of short range and the negativity of the 
pairing field is an indication of a reduced hkelihood of finding two particles at the same location. 

The behaviour of the first, second and third order correlation function is presented in figures 
[TT| [T2] and [13] A generic feature of all three correlation functions is that they become more 
pronounced for smaller particle numbers. For the first order correlation function the diagonal 
has to be identical to one and the deviation in the off-diagonal is fairly small as expected for 
a coherent system. However, the second order density-density correlation is a more sensitive 
probe as this correlation function is less than one, thus exhibits non-classical behaviour. This 
anti-bunching is particularly strong for smaller particle numbers when we approach the Tonks- 
Girardeau regime of a fermionized Bose gas and the correlation function vanishes eventually. 
Recently this effect has been investigated in a number of experiments, e.g. ||5l[6l[T4l and confirms 
the theoretical predictions. The same statements apply to the third order correlation function and 
it can be observed that the deviation from one is even more pronounced. This also implies that 
the third order correlation function ||3l [30l is the most sensitive probe for quantum aspects of the 
field. In addition we notice values which are clearly below one for |y| ^ 1 and x ^ y, because in 



(3) (2) 

this case gx,y ~ gy,y. This can easily be seen by looking at (52 53 1 and taking into account that all 
terms with off-diagonal contributions of the fluctuations in gl^l are negligible for x ^ y. 
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FIG. 12: Second order correlation function gx^y versus x and y. 
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FIG. 14: Correlation functions versus the correlation parameter 7. In subplot a) we depict q and compare 
simulation results for the trapped case (circles) with analytic calculations obtained with Lieb-Liniger theory 

(3) 

(solid line). In subplot b) we depict q and again compare simulation results for the trapped case (circles) 
with analytic calculations obtained with Lieb-Liniger theory (solid line). In both comparisons the circles 
originate for particle numbers N ranging from = 10^ on the left hand side to = 10*^ on the right hand 
side. 

B. Behaviour in the center of the trap 



In figure[T4]we compare the results of our simulations for the second and third order correlation 
function with Lieb-Liniger theory. In contrast to the comparison in subsection IV C| an external 
potential is now included in the calculations with the extended mean-field theory whereas the 
theoretical curve is for a homogeneous gas of bosons. Our simulations are for particle numbers 
ranging from = 10" — 10^ and we only used the values of the correlation functions in the center 
of the trap for the comparison. Compared to subsection IV C| the results in the trapped case deviate 
slightly more from the exact results originating from the homogeneous Lieb-Liniger theory but the 
qualitative behaviour is very similar. 
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FIG. 15: The diagonal second order correlation function gx,x versus 7 for various particle numbers. With 
an increasing value of the correlation parameter 7 the circles correspond to points further outwards from 
the origin. We compare results for the trapped case (circles) with analytic calculations for the second order 
correlation function obtained with LDA-Lieb-Liniger theory (solid line). Plots a) to c) are for N = 10° , 
N = 10^ and N = 10^ (from left to right) and plots d) to f) are for N = 10^, N = 10^ and N = 10^ (from 
left to right). 

C. Diagonal behaviour in the local density approximation 

The local density approximation (LDA) is a frequently employed approximation scheme to 
transfer results of homogeneous systems to spatially trapped gases. It is assumed that a smooth 
variation of the density profile can be incorporated by an adiabatic adjustment of a locally uniform 
gas. The LDA uses a local effective chemical potential [|38l 

fi{x) = fiQ- V{x) = /io - ^mw^x^ , (55) 

where /xq denotes the global equilibrium chemical potential. In order for the LDA to be applicable, 
it is thus necessary that the short-range correlation length is much smaller than the characteristic 
inhomogeneity length. 
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FIG. 16: The diagonal third order correlation function gx,x versus 7 for various particle numbers. With 
an increasing value of the correlation parameter 7 the circles correspond to points further outwards from 
the origin. We compare results for the trapped case (circles) with analytic calculations for the third order 
correlation function obtained with LDA-Lieb-Liniger theory (soUd line). Plots a) to c) are for A?^ = 10° , 

N = 10^ and = 10^ (from left to right) and plots d) to f) are for N = 10^, N = 10^ and N = 10^ (from 
left to right). 

In this context, we want to compare the diagonal behaviour of our numerically calculated cor- 
relation functions to theoretical predictions. By definition, the first order correlation function is 
identical to one along the diagonal and our data behaves accordingly. For the second and third 
order correlation function we will compare our results with the predictions from Lieb-Liniger the- 
ory in the LDA. Naturally the LDA works best in the center of the trap. It can not be expected 
to work in regions where the density drops rapidly and the inhomogeneity length is very small in 
these regions. 

In the Gross-Pitaevskii regime the chemical potential jj, connects the density n to the correlation 
parameter 7, via 

lj,{x) = gn{x), j{x) = , (56) 
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In our simulations we tune the particle number in the trap, which decreases 7 for an increasing 
number of particles. Qualitatively one can expect that the inhomogeneous correlation functions are 
higher than the homogeneous results because in the LDA the external potential leads to a smaller 



chemical potential and according to (56) also to a smaller density compared to the homogeneous 
case. Due to the monotonous decrease of the correlation functions, there is a tendency of the in- 
homogeneous values to be shifted to larger 7 values. All the features that have just been described 



can be seen in figures 15 and 16 where we plotted the correlation functions for particle numbers 
ranging from = 10*^ — 10^ and restricted the plotted regions to the Thomas-Fermi radius. 
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FIG. 17: Finite temperature first order correlation function glc^l versus x and y, for N 
10 huj/kB. 



100 and T 



D. The finite temperature result for a trapped gas 

The zero-temperature results of the previous section can be extended easily to account for finite 
temperature effects [|36ll49ll . One obtains an equilibrium solution for the density matrix G of the 



thermal system (21) from the eigenstates of the self energy matrix (25), according to the Bose- 
Einstein distribution. We present results in the present section for a particle number of = 100 
and a temperature T = 10 huj/kB- 

The main thermal effect is a strong increase of the fluctuations at the edge of the trap at the 
cost of a reduction of the condensate density [|57l . This effect is clearly seen by comparing the 



first order correlation function in figure 17 to the zero temperature result in figure 11 At finite 
temperatures, we also obtain a reduction of first order coherence. Consequently, this leads to 
a situation where the gas is almost thermalized at the edge of the trap, whereas it is coherent 
in the center. The suppression of density fluctuations, also known as anti-bunching, is also less 
pronounced at finite temperature. This can be seen by comparing figures [TS] and [79] to figures [12 



and 13 , which give the zero temperature results. 
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FIG. 18: Finite temperature second order correlation function gx,y versus x and y, for N = 100 and 

T = 10 hio/kB. 




FIG. 19: Finite temperature third order correlation function gx,y versus x and y, for = 100 and T = 

10 hu/kB. 
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FIG. 20: Off-diagonal first 5^^!^. (a) and second order correlation function g^^^_^ (b) versus x. The individual 
curves correspond to temperatures ksT = Ohio (smallest value for = 0) to 10 /lo; (largest value for x = 0) 
with increments of 2 ^w. 



For a thermal gas of noninteracting bosons, one finds gl^l = 2! and Qq^q =3!. It can be seen 



,(3) 



that these values are attained at the edge of the trap where fluctuations dominate. In figure 19 

(3) (2) 

also notice a value of two for \y\ ^ 1 and x y, because we again have gx,y ~ gy,y 



we 



2 in this 



case. 



I I (1 2) 

In figure 20, we present the off-diagonal of the first and second-order correlation function 1^ 
versus x for temperatures from fc^T = — lOhiu with increments of 2hLU. It can be noticed 
that correlations are strongly attenuated with increasing temperature. Looking at the off-diagonal 
we see a reduction of the anti-bunching dip in the center with increasing 



of gS-x in figure 



20 



temperatures, however it is still present for high values of the temperature. It can be understood 
qualitatively from the stronger increase of fluctuations with the temperature at the edge of the trap. 
Thus, the anti-bunching dip in the center of the trap remains visible even at finite temperatures. 
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VI. CONCLUSIONS AND OUTLOOK 

We have presented a detailed study of an extended mean-field theory which shows that it is a 
feasible approach for the description of a weakly correlated gas of bosons in a quasi- ID setup. 
This approach agrees well with exact predictions of zero-temperature Lieb-Liniger theory and can 
easily be applied to spatially inhomogeneous systems at finite temperature. We have not analyzed 
the finite temperature theory of Yang-Yang due its complexity. 

There are many relevant applications for using this extended mean-field theory in different geo- 
metrical configurations like a double- well potential [|58l[59l . Yet another extension of our approach 
is the dimensional crossover out-of-equilibrium where in general the increase of available phase 
space volume leads to a decrease of correlations. An evaluation of the such correlation functions 
is work in progress. 
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APPENDIX A: HIGHER TRANSCENDENTAL FUNCTIONS 

1. Complete elliptic integrals 

Following the definitions and the notation in [|54l . the complete elliptic integral of the first kind 
reads 

/•7r/2 jn 

K{m) = / - (Al) 
"'0 V 1 — m sin 9 



where the parameter < m < 1. For the calculation of mo in section IV A we encounter an 
integral of the form 

where Uc > u). We can show that the evaluation of this integral leads to the complete elliptic 
integral of the first kind by making the substitution k = cot and using m = {ujc — uj)/uJc- 
Similarly the complete elliptic integral of the second kind is defined as 

7r/2 



E{m) = / Vl-msin^^d^. (A3) 



In order to calculate /o in section IV A we end up with the integral 



after separating the constant contribution which leads to the previously discussed integral. The 
same substitution as above k = kc cot simplifies the integral to the form 



h = kj ^ = dd = -K {E{m) - Kim)) . (A5) 

-'0 sin ^^vl — '^sin 6 



2. The Lambert-W function 

The Lambert- VT-function is implicitly defined by the solution of the transcendental equation 

m 

z = We^. (A6) 
In the case of large arguments z ^ 1, one can use an asymptotic expansion 

w{z) = ei-e2 + i2/ii + ..., (A7) 

with ii = In 2; and £2 = In In z. 
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FIG. 21: Integration contour for the evaluation of (Bl i. 



APPENDIX B: DEFORMATION OF THE INTEGRATION CONTOUR IN THE COMPLEX 

PLANE 



The results of subsection IV B 2 can be derived alternatively with the help of complex inte- 
gration [[60ll . If we take the inverse Fourier transform of the anomalous fluctuation of (38), we 
get 



-ikr 



dk . 



Using the substitutions k = kz, r' = kr and = k'l/k'^ this equation reduces to 

e~*^'^ dz 



(Bl) 



m[r] 



(B2) 



, Vz^ + b^^MTT k 

For the evaluation of this integral we make a branch cut between —i and —ib and choose the path 



of integration as can be seen in Fig. 21 



The contributions from Ci and Cq vanish if the contour is moved to infinity and the contribu- 
tions from C2 and C5 cancel each other. The integrals along the semicircles around —ib and the 
circle around —i tend to zero if the radius tends to zero. Due to the branch cut the contributions 
from C3 and C4 are equal. Thus, the integral to be solved reads 

J _ e-'^-'^dz _ r' 2e-^'-'^dz 



(B3) 
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and by changing the variable of integration (z — —i — iy), taking into account that 6 » 1 in the 
Gross-Pitaevskii regime, we get 

I ^ 2e-' r , (B4) 
Jo - 1 - 27/ - y2 

As we are looking for an approximation for r' = ^ 1, we notice that only small values of y 
play an important role for the evaluation of the integral. 

Hence we neglect the expression 1 + 2y + y'^ va the second term in the denominator which 

yields 

-r' roo -yr' 9 

and as m is an even function in r we get the final result 

m(r) ^ -^K^{k\r\) « Mk\r\) . (B6) 
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